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ABSTRACT 

A  previously  derived  approximation  to  the  standard  Fourier  integral  technique  for  linear  mountain  waves  is 
extended  to  include  nonhydrostatic  effects  in  a  background  flow  with  height-dependent  wind  and  stratification. 
The  approximation  involves  using  ray  theory  to  simplify  the  vertical  eigenfunctions.  The  generalization  to 
nonhydrostatic  waves  requires  special  treatment  for  resonant  modes  and  caustics.  Resonant  modes  are  handled 
with  a  small  amount  of  damping,  and  caustics  are  handled  with  a  uniformly  valid  approximation  involving  the 
Airy  function.  This  method  is  developed  for  both  two-  and  three-dimensional  flows,  and  its  results  are  shown 
to  compare  well  with  an  exact  analytical  result  for  two-dimensional  mountain  waves  and  with  a  numerical 
simulation  for  two-  and  three-dimensional  mountain  waves. 


1.  Introduction 

The  Fourier  integral  representation  has  been  used  in 
the  linear  modeling  of  mountain  waves  for  over  50  yr 
and  remains  an  important  tool  for  mountain  wave  theory 
and  forecasting  (e.g.,  Vosper  and  Mobbs  1996;  Smith 
et  al.  2002;  Shutts  1997).  The  representation  is  ex¬ 
pressed  as  the  inverse  Fourier  transform,  over  the  hor¬ 
izontal  wavenumbers,  of  the  vertical  eigenfunctions.  Ex¬ 
cept  for  very  special  cases,  the  vertical  eigenfunctions 
must  be  calculated  numerically  as  the  solution  of  a  two- 
point  boundary  value  problem.  The  numerical  procedure 
is  complicated  by  the  presence  of  critical  layers  and 
trapped  waves  (e.g.,  Simard  and  Peltier  1982;  Shutts 
1997). 

In  Broutman  et  al.  (2002)  the  authors  simplified  the 
calculation  by  approximating  the  vertical  eigenfunctions 
with  ray  theory.  We  considered  a  slowly  varying  height- 
dependent  background  in  the  hydrostatic  limit.  Here,  we 
extend  this  work  to  the  nonhydrostatic  case,  which  re¬ 
quires  the  treatment  of  trapped  waves  and  the  associated 
ray  singularity  at  caustics. 

The  ray  approximation  for  the  vertical  eigenfunctions 
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is  not  to  be  confused  with  the  ray  approximation  for  the 
spatial  wavefield.  The  latter  has  been  computed  for 
mountain  waves  by  Broad  (1999)  and  Shutts  (1998), 
among  others,  and  can  be  obtained  in  principal  by  the 
method  of  stationary  phase.  This  is  equivalent  to  taking 
an  inverse  Fourier  transform  of  the  vertical  eigenfunc¬ 
tions  and  then  finding  the  ray  approximation.  Here,  we 
switch  the  order  of  these  operations.  That  is,  we  find 
the  ray  approximation  for  the  vertical  eigenfunctions 
and  then  take  an  inverse  Fourier  transform.  The  order 
of  operations  makes  an  important  difference  for  the  so¬ 
lution  near  caustics. 

A  caustic  is  a  singularity  of  ray  theory  that  occurs 
where  neighboring  rays  intersect  each  other,  an  extreme 
violation  of  the  slowly  varying  approximation.  For 
trapped  waves,  there  can  be  numerous  caustics  in  the 
lee  of  the  mountain,  which  is  why  we  do  not  calculate 
the  spatial  ray  solution  in  the  present  paper.  While  these 
caustics  are  correctable  with  an  Airy  function,  typically, 
the  correction  requires  quantities  that  are  not  easily  ob¬ 
tained  unless  the  ray  solution  is  known  analytically. 

The  ray  approximation  for  the  vertical  eigenfunctions 
also  has  caustics,  and  these  too  are  correctable  with  an 
Airy  function,  typically.  But  the  caustics  for  the  vertical 
eigenfunctions  and  the  caustics  for  the  spatial  ray  ap¬ 
proximation  have  different  properties.  For  example,  the 
caustics  for  the  vertical  eigenfunctions  always  occur  at 
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Fig.  1.  Ray  paths  for  Long’s  model  all  lie  on  the  vertical  axis  above  the  mountain  when 
plotted  in  the  spatial  domain  (left)  but  are  separated  by  horizontal  wavenumber  k  when 
plotted  in  the  Fourier  domain  (right). 


a  buoyancy  frequency  turning  point,  while  the  caustics 
for  the  spatial  ray  approximation  never  occur  at  a  buoy¬ 
ancy  frequency  turning  point.  We  explain  why  and  show 
an  example  in  section  3 — see  the  description  of  Fig.  2. 

Most  importantly,  the  correction  of  the  caustics  for 
the  vertical  eigenfunctions  is  a  one-dimensional  anal¬ 
ysis.  This  simplifies  the  calculation,  compared  with  the 
multidimensional  case  for  the  spatial  caustics.  Thus,  the 
procedure  we  use  is  to  compute  the  ray  approximation 
for  the  vertical  eigenfunctions,  correcting  for  caustics 
in  a  standard  way  with  the  Airy  function.  An  inverse 
Fourier  transform  of  the  result  then  gives  the  spatial 
solution,  which  is  valid  at  all  of  the  spatial  caustics 
without  further  correction. 

The  present  approach  is  related  to  Maslov’s  method, 
developed  in  the  1960s  by  V.  P.  Maslov  (e.g.,  Maslov 
and  Fedoriuk  1981;  Kravtsov  and  Orlov  1999).  In  Mas¬ 
lov’s  method,  caustics  in  the  spatial  domain  are  avoided 
by  using  the  inverse  Fourier  transform  of  the  ray  so¬ 
lution  from  the  Fourier  domain.  In  our  problem,  the  ray 
solution  from  the  Fourier  domain  is  just  the  ray  solution 
for  the  vertical  eigenfunctions.  The  hydrostatic  moun¬ 
tain  wave  theory  of  Shutts  (1998)  is  an  example  of  this 
approach,  as  we  discussed  in  Broutman  et  al.  (2002). 
We  have  since  noticed  that  a  similar  idea  was  suggested 
for  mountain  waves  by  Miles  (1969),  though  the  rele¬ 
vant  expression  given  by  Miles  in  his  Eq.  (4.13)  does 
not  apply  to  trapped  waves. 

The  strength  of  Maslov’s  method  is  how  it  deals  with 
nonseparable  problems,  which  lack  an  eigenfunction 
representation.  Maslov’s  method  combines  ray  solutions 
from  the  Fourier  and  spatial  domains  using  Fourier 
transforms,  in  such  a  way  that  the  local  caustic  correc¬ 
tion  functions,  like  the  Airy  function,  are  never  needed. 
Brown  (2000)  gives  an  example.  The  same  procedure 
could  be  used  here,  but  since  our  problem  has  an  ei¬ 
genfunction  representation,  it  is  natural  to  use  an  Airy 
function  for  the  caustics  in  the  Fourier  domain  because, 
as  we  noted  above,  the  analysis  in  the  Fourier  domain 
is  a  simple  one-dimensional  calculation. 

One  way  to  derive  the  ray  approximation  for  the  ver¬ 
tical  eigenfunctions  is  by  an  asymptotic  analysis  of  the 
vertical  structure  equation,  given  in  general  form  and 
in  several  limits  by  Inverarity  and  Shutts  (2000).  This 
is  the  approach  used  by  Shutts  (1998)  and  Miles  (1969). 
Another  derivation,  which  is  used  here  and  leads  to  the 


same  result,  is  to  ray  trace  in  the  Fourier  domain,  as  we 
did  in  Broutman  et  al.  (2002)  for  the  hydrostatic  case. 
This  approach  is  applicable  to  more  general  back¬ 
grounds,  for  example,  with  horizontal  variability,  but 
mostly  we  use  it  because  such  an  alternate  derivation 
helps  to  enhance  the  understanding  of  the  results. 

To  analyze  trapped  waves,  it  is  convenient  to  start 
with  the  case  of  uniform  stratified  hydrostatic  flow  over 
topography  in  a  channel  (section  2).  The  trapping  is 
imposed  by  the  rigid  upper  lid.  This  is  Long’s  model, 
and  the  present  analysis  leads  to  a  new  derivation  of 
the  linearized  limit  of  Long’s  solution.  We  can  easily 
extend  the  derivation  to  the  case  of  nonhydrostatic 
waves  in  a  height-dependent  background  without  an  up¬ 
per  boundary.  This  is  done  in  section  3. 

In  section  4  we  show  that  the  present  method  com¬ 
pares  well  with  the  two-dimensional  theory  of  Wurtele 
et  al.  (1987)  based  on  the  exact  vertical  eigenfunctions. 
In  section  5  we  give  a  three-dimensional  example. 
Throughout  the  paper  we  use  the  Boussinesq  approxi¬ 
mation,  and  we  ignore  the  effects  of  the  earth’s  rotation. 

2.  Long’s  hydrostatic  channel  model 

A  simple  demonstration  of  the  present  method  for  the 
case  of  trapped  waves  is  provided  by  Long’s  model  of 
flow  over  topography  in  a  channel.  Long’s  model  (Long 
1955;  Baines  1995)  is  hydrostatic  and  steady  state,  with 
constant  mean  flow  U  and  constant  mean-buoyancy  fre¬ 
quency  N.  The  topography  z  =  h(x)  is  localized,  with 
h  — >  0  as  x  ±oo.  The  upper  boundary  of  the  channel 
is  at  z  =  D. 

a.  Rays  in  the  spatial  domain 

The  rays  in  the  spatial  domain  all  emerge  from  a 
single  point  at  the  center  of  the  topography.  [This  point 
source  of  rays  is  a  feature  of  radiation  from  a  localized 
source  and  can  be  justified  by  a  stationary-phase  ap¬ 
proximation,  as  in  Shutts  (1998).]  The  rays  then  prop¬ 
agate  straight  upward  along  the  vertical  axis  above  the 
mountain  (Lig.  1,  left).  The  rays  cannot  spread  laterally 
in  the  hydrostatic  limit  because  the  upstream  intrinsic 
group  speed  of  the  waves  is  perfectly  balanced  by  the 
downstream  flow.  Nor  can  the  rays  reflect  laterally  from 
the  mountain  after  a  round  trip  to  the  top  of  the  channel 
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and  back.  Even  if  the  mountain  were  sloped  at  the  point 
where  the  rays  return  to  it,  the  rule  for  internal  waves 
is  that  the  rays  preserve  their  angle  with  respect  to  the 
vertical  upon  reflection.  Thus  the  rays  reflect  straight 
upward  again.  Without  any  lateral  dispersion,  the  slowly 
varying  approximation  fails  completely  and  the  spatial 
ray  solution  predicts  a  vertical  displacement  amplitude 
that  is  infinite  at  any  point  on  the  vertical  axis  and  that 
is  zero  at  any  point  off  the  vertical  axis. 

b.  Rays  in  the  Fourier  domain 

We  can  make  more  progress  with  ray  theory  if  we 
express  the  ray  solution  as  a  function  of  the  horizontal 
wavenumber  k  and  the  depth  z.  This  is  the  Fourier  do¬ 
main  for  this  problem.  The  depth  is  treated  parametri¬ 
cally  and  is  not  transformed. 

For  Fong’s  model,  the  rays  that  line  up  on  the  vertical 
axis  above  the  mountain  in  the  spatial  domain  are  sep¬ 
arated  in  the  Fourier  domain  by  their  values  of  k,  as 
illustrated  in  the  right  side  of  Fig.  1.  There  is  still  an 
infinite  number  of  rays  of  the  same  k  at  each  z,  due  to 
the  continual  generation  of  rays  at  the  lower  boundary 
and  the  continual  reflection  of  rays  by  both  boundaries. 
It  is,  however,  possible  to  sum  the  contributions  from 
all  rays  of  the  same  k,  as  we  now  show. 

Fet  u>  be  the  intrinsic  frequency,  with  a>  =  — kU  for 
stationary  waves.  We  assume  that  U  is  positive  and  that 
k  is  negative.  The  hydrostatic  dispersion  relation  gives 
the  vertical  wavenumber  magnitudes  as  m  =  \k\N/w, 
which  after  substituting  for  a>  becomes  m  =  N/U. 

The  ray  solution  for  the  vertical  displacement  in  the 
Fourier  domain  will  be  represented  by  fj(k,  z).  The 
boundary  conditions  are 

fj(k,  0)  =  h(k),  t )(k,  D)  =  0,  (1) 

where  h(k )  is  the  Fourier  transform  of  h(x). 

We  derive  the  ray  solution  for  fj(k,  z)  by  tracing  rays 
in  the  Fourier  domain.  We  start  with  some  initial  ray 
generated  at  z  =  0  and  ignore  all  rays  generated  previ¬ 
ously.  That  is,  we  consider  an  initial  value  problem  with¬ 
out  the  complication  of  transients  and  look  for  the  long¬ 
time  solution.  Just  before  the  initial  ray  reaches  the  upper 
boundary  for  the  first  time,  the  ray  solution  is 


Just  before  the  initial  ray  reaches  the  upper  boundary 
for  the  second  time,  after  one  reflection  from  the  upper 
boundary  and  one  reflection  from  the  lower  boundary, 
the  ray  solution  is 

fj(2>  =  h(k)e~imz  —  h(k)eimU~2D)  +  h(k)e~im(z+2D\  (4) 

The  third  term  is  the  same  as  the  first  term  except  for 
the  phase  shift  of  —m2D,  resulting  from  a  difference  in 
the  length  of  the  propagation  path. 

Each  incident  ray  combines  with  its  reflected  ray  to 
nullify  the  vertical  displacement  at  the  boundary.  The 
upper  boundary  condition  is  thus  satisfied  pairwise.  The 
lower  boundary  condition  is  thus  satisfied  by  the  single 
ray  being  generated  for  the  first  time  at  the  mountain 
and  represented  in  the  first  term  on  the  right  in  (2),  (3), 
and  (4). 

The  continued  reflections  generate  an  infinite  number 
of  rays  that  superimpose  at  each  z .  We  represent  this 
for  the  upward-moving  rays  by 

z)  =  he-‘mzS,  (5) 

and  for  the  downward-moving  rays  by 

T7dn (k,  z)  =  —heim(z^2D)S,  (6) 

with 

5  =  2  e~in2mD  (7) 

n=0 

For  a  fixed  z,  each  term  in  the  sum  for  rjup  and  each 
term  in  the  sum  for  fjdn  differs  only  by  a  phase  shift, 
an  integer  multiple  of  —m2D. 

The  sum  5  is  divergent  in  the  strict  sense,  since  the 
terms  do  not  vanish  as  n  — >  °o,  but  it  does  converge  in 
the  sense  of  generalized  functions.  From  Hardy  [1949, 
their  Eq.  (1.2.2)], 

5  =  (1  —  i  cotmD)/2  =  —  ieimD/2  sin mD.  (8) 

The  total  vertical  displacement  thus  becomes 
fj(k,  z)  =  rjup  +  fjdn 

=  hg-imD  Yg  —  im (z~D)  —  gim(z—D) 

_  h  sin m(D  -  z) 
sin  mD 


fjm  =  h(k)e-‘mz.  (2) 

The  zero  superscript  indicates  that  no  reflections  have 
occurred  yet.  With  notation  m  >  0,  the  negative  sign 
in  the  phase  indicates  a  downward  phase  velocity  and 
an  upward  group  velocity.  After  the  initial  ray  has  its 
first  reflection  at  the  top  boundary,  but  just  before  it  has 
its  second  reflection  at  the  bottom  boundary  the  ray 
solutions  is 


The  spatial  solution  is  now  obtained  by  inverse  Fourier 
transform: 


Vix,  z) 


7)(k,  z)eikx  dk 


sin  m(D  —  z) 
sin  mD 


h(k)eikx  dk 


fjm  =  h(k)e  imz  —  h(k)eim<~z  2D).  (3) 

The  second  term  represents  downward-moving  rays  and 
combines  with  the  first  term  to  satisfy  the  upper  bound¬ 
ary  condition. 


,  ,  N  sinm(D  -  z) 

=  li(x) - 

sinwiD 


(10) 


This  is  Fong’s  solution  for  small-amplitude  topography 
(see  p.  279  of  Baines  1995). 
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3.  Two-dimensional  theory  for  U(z),  N(z) 

We  now  remove  the  rigid  lid  of  section  2  and  examine 
a  two-dimensional  model  in  which  the  waves  are  non¬ 
hydrostatic  and  trapped  by  variations  in  the  background 
wind  and  stratification.  The  generalization  to  three  di¬ 
mensions  is  straightforward  (section  5).  We  develop  the 
theory  for  a  height-dependent  mean  wind  U{z)  and  mean 
buoyancy  frequency  N{z).  The  nonhydrostatic  disper¬ 
sion  relation  is 


i b  =  ±kN/(k2  +  m2)m.  (11) 

For  stationary  waves  &  =  —kU.  We  assume  slowly 
varying  waves,  except  near  caustics.  This  limits  the 
strength  of  variations  in  U(z),  N(z)  (see,  e.g.,  section  3c 
of  Broutman  et  al.  2002). 

As  in  section  2,  we  use  an  infinite  sum  to  represent 
the  multiply  reflected  rays  in  the  Fourier  domain.  The 
theory  of  section  2  needs  to  be  modified  for  the  present 
model  to  account  for  refraction  in  a  height-dependent 
background  and  to  account  for  the  phase  shift  from  a 
buoyancy  frequency  turning  point. 

First  we  consider  the  wave  amplitude.  In  section  2, 
each  ray  has  a  constant  vertical-displacement  amplitude 
h(k)  in  the  Fourier  domain — see  Eqs.  (5)-(7).  For  a 
height-dependent  background,  this  must  be  modified  to 
conserve  the  vertical  flux  of  wave  action  for  each  k. 
Defining 


G  =  N2cgJu>,  (12) 

the  wave-action  flux  is  Gif) I2.  Enforcing  the  lower 
boundary  condition  f)  =  h  at  z  =  0  gives 

If)  I  =  /z[G0/G]1/2,  (13) 

where  G0  =  G(z  =  0). 

Next  we  consider  the  wave  phase.  The  channel  height 
D  in  Long’s  model  is  replaced  by  the  turning  point 
height  z„  which  is  a  function  of  k.  The  generalization 
of  the  wave  phase  in  (5)-(6)  is  the  standard  one  in  which 
mz  becomes  the  phase  integral  \'0m  dz,  and  accordingly 
mD  becomes  jj )m  dz.  There  is  a  phase  shift  at  the  turning 
point  z„  and  this  enters  into  the  generalization  of  S  in 
(7),  as  we  now  discuss. 

Following  a  ray  through  one  round-trip  to  the  top  of 
the  channel  and  back  led,  in  section  2,  to  a  decrease  in 
the  wave  phase  by  an  amount  —2 mD.  Actually  there 
are  phase  shifts  of  77  at  both  the  upper  and  lower  bound¬ 
aries,  but  this  contributes  an  additional  factor  of  277  to 
the  phase  of  S  and  so  can  be  ignored.  In  the  present 
generalization  we  still  have  the  phase  shift  of  77  at  the 
lower  boundary,  but  the  upper  boundary  is  a  buoyancy- 
frequency  turning  point,  and  this  introduces  a  different 
phase  shift  which  we  call  a,  for  now.  Thus  2 mD  in  (7) 
becomes  2$o'm  dz  +  a  +  77. 

Summarizing,  the  generalization  of  the  ray  solution 
from  Long’s  model  is  then,  in  (5)-(6), 


h  -A  /z[G0/G]1/2, 


mz  -A  m  dz, 


mD  -a  m  dz, 


and  in  (7), 


2mD  -A  2  m  dz  +  a  +  tt. 


(14) 

(15) 

(16) 


(17) 


Making  these  substitutions,  and  continuing  in  the 
same  way  as  the  derivation  that  leads  to  (9),  gives 


where 


z) 


sin[r  -  x(z)] 
sinT 


(18) 


X 


m(zf)  dz1, 


(19) 


r  =  X(zt)  +  a/2  +  77/2.  (20) 


Long’s  solution  for  f)  in  (9)  is  recovered  from  (18) 
by  setting  a  =  77  (the  phase  shift  for  reflection  from  a 
solid  boundary),  and  by  letting  x  =  mz  an(i  G/Gn  =  1 
(as  they  are  for  a  uniform  medium). 

Equation  (18)  applies  only  to  those  k  values  that  have 
turning  points.  Lor  the  k  values  without  turning  points 
we  use  the  formula  derived  in  Broutman  et  al.  [2002, 
their  Eq.  (26)].  In  present  notation,  this  is 


7j(k,  z)  =  h(G0/G)ln  exp[-z 


m(k,  z’)  dz’],  (21) 


which  represents  only  upward-moving  rays. 

The  inverse  Fourier  transform  (10)  of  (18),  and  pos¬ 
sibly  of  (21),  gives  the  spatial  solution  for  p(x,  z).  But 
there  are  two  singularities  to  deal  with  before  this  can 
provide  a  meaningful  result.  The  first  is  G  =  0,  which 
occurs  when  cg3  =  0  and  corresponds  to  either  a  critical 
layer  or  a  turning  point.  We  dealt  with  the  critical  layer 
in  Broutman  et  al.  (2002)  by  simply  filtering  out  Fourier 
components  of  large  vertical  wavenumber.  That  worked 
effectively,  at  least  for  the  model  considered  in  that 
paper.  Now  we  must  deal  with  the  turning  point  sin¬ 
gularity,  which  is  a  caustic  in  the  Fourier  domain.  The 
other  singularity  of  (18)  is  F  =  nir,  for  integer  n,  which 
occurs  when  reflecting  rays  of  the  same  k  have  perfect 
constructive  interference. 


a.  Caustics 

The  standard  correction  procedure  for  a  caustic  is  to 
match  the  ray  solution  to  an  Airy  function  Ai(p)  in  its 
asymptotic  form  for  large  negative  p  (Bender  and  Or- 
szag  1978): 
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Ai(p) 


2\pl 


2  77 

_|p|3/2  +  _ 

3  4 


(22) 


The  caustic  height  z  —  z,  corresponds  to  p  =  0,  with  p 
>  0  (z  >  z,,)  in  the  “shadow”  zone  and  p  <  0  (z  <  zt) 
in  the  “lit”  zone.  The  relation  between  p  and  z  is 


2 

-lpl3/2  = 


=  r 


m(z')  dz!  =  x(z,)  -  X(z) 

1 

-  -(a  +  77)  -  *(z). 


(23) 


Setting  the  previously  undetermined  phase  shift  at  the 
caustic  to  a  =  7r/2  gives 


|ipi3/2  +  j  =  r  -  *(z), 


(24) 


so  the  sine  functions  in  the  Airy  function  formula  (22) 
and  in  the  ray  solution  (18)  have  the  same  argument, 
as  required. 

The  matching  of  the  ray  solution  to  the  Airy  function 
is  a  well-known  procedure,  but  we  will  use  an  alternate 
description  that  gives  a  uniform  approximation,  that  is, 
a  single  expression  valid  close  to  and  far  from  the  caus¬ 
tic.  This  approach  is  generally  more  accurate  than  the 
matching  method.  A  discussion  of  the  uniform  approx¬ 
imation  for  the  one-dimensional  Helmholtz  equation  is 
given  in  Bender  and  Orzsag  (1978,  p.  510),  and  for  a 
more  general  case,  by  Kravtsov  and  Orlov  (1999).  The 
uniform  solution,  derived  in  the  appendix,  is 


Vu  =  ^ 


Glp 


G2p, 


Ai  (p) 
Ai(po) 


(25) 


The  zero  subscript  means  an  evaluation  at  z  =  0.  The 
factor  of  G~m  causes  the  ray  solution  (18)  to  be  singular 
at  the  caustic,  but  in  (25)  that  singularity  is  canceled  by 
the  factor  I  pi 1/4 .  The  spatial  solution  is  now 


V(x,  z) 


fj„(k,  z)eihc  dk.  (26) 


b.  Resonant  modes,  F  =  inr 

In  the  Fourier  domain,  the  total  phase  change  2 T 
along  the  ray  in  one  round-trip  from  the  ground  to  the 
turning  point  and  back  again  (including  the  phase  shifts 
at  the  caustic  and  at  the  ground)  can  turn  out  to  be  an 
integer  multiple  of  2tt.  We  then  have  a  resonant  mode 
in  which  there  is  perfect  constructive  interference  be¬ 
tween  all  rays  of  the  same  k.  In  Long’s  model,  all  modes 
have  the  same  resonant  condition  mD  =  n  it  for  integer 
n — see  (10).  There  is  no  linear  steady-state  solution  in 
that  case.  In  the  present  model,  the  resonant  condition 
is  different  for  each  k,  and  each  gives  rise  to  an  inte¬ 
grate  pole  in  the  inverse  Fourier  transform,  as  in  Eq. 
(7)  of  Wurtele  et  al.  (1989). 


We  remove  the  singularity  associated  with  the  reso¬ 
nant  modes  by  adding  a  small  amount  of  damping  in 
the  form  of  an  imaginary  part  for  k.  We  discuss  the 
effectiveness  of  this  damping  in  the  next  section.  (See 
the  description  of  Fig.  4.) 

4.  Comparison  with  Wurtele  et  al.  (1987) 

Wurtele  et  al.  (1987),  hereafter  WSK,  have  analyzed 
two-dimensional  mountain  waves  in  a  mean  wind  given 
by 

U  =  U0Z,  Z  =  1  +  z/L,  (27) 

where  U0  and  L  are  constants.  We  consider  the  tropo¬ 
spheric  model  in  WSK  section  2a.  The  topography  used 
by  WSK  is 

h(x)  =  h„a2/(x2  +  a2),  (28) 

where  a  is  a  constant.  The  Fourier  transform  of  h  is 
h  =  h0ae-"lkl  (29) 

In  the  following  we  use  the  parameter  values  of 
WSK’s  case  II,  which  are  N  =  0.01  s_1  and 

a  =  2.5  km,  h0  =  100  m,  L  =  4.0  km, 

U0  =  10  m  s_1.  (30) 

Figure  2  shows  ray  paths  for  this  model,  plotted  as 
a  function  of  their  spatial-domain  coordinates  (top)  and 
as  a  function  of  their  Fourier-domain  coordinates  (bot¬ 
tom).  The  spatial  rays  originate  at  the  origin  and  are 
calculated  from  an  integration  of  the  standard  ray  equa¬ 
tions  (Lighthill  1978).  The  corresponding  rays  in  the 
two  domains  can  be  identified  by  the  height  of  their 
turning  points,  which  are  the  same  in  the  two  domains. 
In  the  Fourier  domain  the  rays  follow  a  straight  vertical 
line  because  k  is  constant  along  the  ray  in  a  horizontally 
uniform  background. 

In  the  Fourier  domain,  the  turning  point  is  a  caustic. 
Neighboring  rays  of  the  same  k,  generated  at  slightly 
different  times,  intersect  as  they  reflect  from  the  turning 
point.  The  intersection  of  neighboring  rays  is  what  de¬ 
fines  a  caustic.  In  the  spatial  domain  the  turning  point 
is  not  a  caustic.  This  is  most  clear  at  the  first  downwind 
occurrence  of  turning  points  in  the  top  of  Fig.  2.  The 
rays  arc  back  toward  the  ground  without  intersecting. 
Caustics  form  for  a  particular  ray  only  after  that  ray 
reflects  once  from  the  ground.  Four  caustics  are  dis¬ 
cernible  in  the  top  of  Fig.  2  as  approximately  straight 
lines  that  slope  upward  with  increasing  x.  Each  caustic 
is  an  envelop  for  the  rays  that  intersect  it,  and  each  ray 
intersects  its  caustic  before  reaching  its  turning  point. 

The  turning  point  in  the  spatial  domain  is  never  a 
caustic  in  this  model.  To  understand  why,  note  that  each 
ray  is  uniquely  determined  by  its  position  ( x ,  z)  and  by 
its  horizontal  wavenumber  k.  At  a  caustic  in  the  spatial 
domain,  the  intersecting  rays  must  by  definition  have 
the  same  value  of  x  and  z.  At  a  caustic  in  the  Fourier 
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Fig.  2.  Corresponding  ray  paths  in  the  spatial  domain  (top)  and  in  the  Fourier  domain 
(bottom)  for  the  model  of  section  4.  The  values  for  the  horizontal  wavenumber  k  correspond 
to  horizontal  wavelengths  ranging  from  6.5  to  30  km.  In  the  bottom,  the  k  axis  is  normalized 
by  the  mountain  width  a  =  2.5  km,  and  the  dotted  line  indicates  the  turning-point  height. 


domain,  the  intersecting  rays  must  by  definition  have 
same  value  of  k  and  z.  If  the  ray  intersections  were  to 
happen  at  the  same  z  in  both  domains,  the  intersecting 
rays  would  then  have  the  same  values  of  x,  z,  and  k, 
and  thus  would  not  be  distinct. 

The  exact  vertical  eigenfunctions  for  the  WSK  model 
involve  Bessel  functions,  as  contained  in  WSK’s  Eq. 
(6).  Converting  their  expression  from  perturbation 
streamfunction  to  vertical  displacement  and  removing 
their  non-Boussinesq  mean-density  scaling,  the  WSK 
eigenfunctions  become 

K,JkLZ ) 

W*.  =  hf,aZ,ne  .  (31) 

KJkL) 

Here,  /r  =  N2L2/U20  —  1/4,  and  Kifl  is  the  modified 
Bessel  function  of  imaginary  order. 

Our  f)u{k,  z)  in  (25)  is  an  approximation  for  the  exact 
vertical  eigenfunction  (31).  Tests  of  this  approximation 
are  shown  in  Fig.  3,  in  which  we  plot  the  exact  eigen¬ 
function  (31),  the  uniform  approximation  (25),  and  the 
singular  ray  solution  (18).  Each  of  the  three  panels  cor¬ 
responds  to  a  different  value  of  k,  as  labeled  by  the 
horizontal  wavelength  A  =  2ir/\  k\ .  The  agreement  be¬ 
tween  the  exact  vertical  eigenfunction  and  the  uniform 
approximation  is  close,  with  the  case  of  A  =  40  km 
showing  the  largest  discrepancy.  The  size  of  the  dis¬ 
crepancy  is  related  to  the  nearness  of  the  vertical  ei¬ 
genfunction  to  a  resonant  mode.  At  resonance  the  ver¬ 
tical  eigenfunction  diverges,  and  the  difference  between 


the  exact  and  approximate  values  is  accentuated,  as  we 
show  next. 

In  Fig.  4,  the  uniform  approximation  and  the  exact 
vertical  eigenfunction  are  plotted  as  a  function  of  \k\a 
at  a  height  corresponding  to  p  =  —1,  close  to  the  depth 
at  which  the  vertical  eigenfunctions  have  their  local 
maximum  nearest  the  caustic.  Resonance  occurs  at  A 
approximately  equal  to  15  and  36  km,  or  \k\a  ~  1.05 
and  \k\a  ~  0.44,  respectively. 

The  uniform  approximation  in  Fig.  4  is  calculated 
with  an  imaginary  part  A,  =  5  X  10~6  m_1  added  to  the 
real  part  of  k.  The  only  significant  effect  of  A,  is  on  those 
eigenfunctions  that  are  near  resonance.  The  asterisks 
show  the  three  values  of  ka  corresponding  to  Fig.  3, 
which  was  also  computed  with  ki  =  5  X  10~6.  From 
left  to  right  in  Fig.  4,  the  asterisks  denote  the  values  A 
=  40,  20,  10  km.  The  accuracy  of  the  uniform  approx¬ 
imation  in  Fig.  4  is  seen  to  depend  on  the  nearness  of 
A  to  a  resonant  value. 

Note  that  by  damping  the  near-resonant  modes  pref¬ 
erentially,  as  indicated  by  Fig.  4,  we  damp  the  waves 
of  the  far  field  preferentially.  In  the  far  field,  the  resonant 
modes  make  the  dominant  contribution  to  the  wave  field, 
as  noted  in  WSK’s  discussion  of  their  Eq.  (7).  In  the 
near  field,  a  wide  range  of  eigenfunctions  contributes 
significantly  to  the  solution.  Thus  damping  the  near¬ 
resonant  modes  has  less  of  an  effect  on  the  near-field 
solution  than  on  the  far- field  solution. 

The  value  of  ki  =  5  X  10~6  m_1  is  used  to  generate 
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Fig.  3.  A  comparison  of  the  exact  vertical  eigenfunction  (31),  the  uniform  approximation  (25), 
and  the  ray  approximation  (18).  The  solutions  are  normalized  by  h(k).  Three  cases  are  given  for 
different  values  of  A  =  2tt/\I<  \  .  The  vertical  dotted  line  indicates  the  turning-point  height,  where 
the  ray  approximation  diverges.  Parameter  values  are  given  in  (30). 


the  results  plotted  in  this  paper  and  corresponds  to  a 
horizontal  decay  scale  ky1  of  200  km.  We  chose  this 
value  specifically  for  the  WSK  model,  in  which  the 
mountain  width  is  only  2.5  km.  More  generally,  we 
would  choose  the  value  of  k in  order  to  damp  the  far- 
held  waves  sufficiently  so  that  periodic  wrap-around 
errors  are  not  important.  Periodic  wrap-around  errors 
result  from  the  approximation  of  the  inverse  Fourier 
transform  (26)  by  a  finite  Fourier  series  (as  described 
in,  e.g.,  Broutman  et  al  2002).  Note  that  periodic  wrap¬ 
around  errors  in  the  discrete  approximation  of  the  in¬ 
verse  Fourier  transform  also  occur  when  the  exact  ver¬ 
tical  eigenfunctions  are  used  (e.g.,  Shutts  2001).  For  the 
two-dimensional  case  in  Fig.  5,  to  be  discussed  next, 
we  found  that  changing  ky'  from  200  to  1000  km  in¬ 
creased  the  maximum  wave  amplitudes  by  about  15%. 


For  the  three-dimensional  case  of  Fig.  6,  the  same 
change  in  ky'  increased  the  maximum  wave  amplitudes 
by  about  5%. 

In  Fig.  5,  we  plot  the  solution  for  the  vertical  velocity 
predicted  by  the  present  simplified  Fourier  method.  This 
is  computed  using  the  ray-approximated  vertical  eigen¬ 
function  w(k,  z)  =  ikUrj„(k,  z).  The  result  agrees  well 
with  the  result  of  a  numerical  model,  which  is  also 
shown  in  Fig.  5.  The  numerical  model  is  a  standard  one: 
it  is  nonlinear,  Boussinesq,  and  discretized  with  second- 
order  finite  differences.  The  lower  boundary  is  free  slip, 
and  the  flow  is  turned  on  instantly  to  full  strength  at 
the  initial  time.  The  results  are  plotted  after  15  h  of 
simulation  time,  when  an  approximate  steady  state  is 
achieved.  The  numerical  model  has  a  domain  size  of 
130  km  by  52  km,  with  a  grid  spacing  of  1  km  in  each 
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Fig.  4.  The  exact  vertical  eigenfunction  (31)  and  the  uniform  approximation  (25)  as  a  function 
of  I  k  I  a  at  a  depth  corresponding  to  the  Airy  function  variable  p  =  —  1 .  The  uniform  calculation 
has  an  imaginary  wavenumber  =  5  X  10-6  m-1.  The  exact  calculation  is  for  =  0.  The 
left,  center,  and  right  asterisks  correspond  respectively  to  A  =  40,  20,  10  km,  the  values  used 
in  Fig.  3. 


direction.  A  sponge  is  in  effect  from  about  x  >  80  km 
and  z  >  35  km.  The  results  from  this  numerical  model 
are  similar  to  those  obtained  by  numerical  simulation 
in  Fig.  9  of  WSK.  [Our  magnitudes  are  lower  than  those 
plotted  in  Figs.  2  and  9  of  WSK  by  a  factor  of  about 
5.  We  have  concluded  that  there  is  a  minor  error  in  the 
labeling  of  these  two  figures  in  WSK  since  we  have 
verified  our  results  with  two  different  numerical  models 
and  an  evaluation  of  WSK’s  far-held  approximation  (7).] 


5.  A  three-dimensional  example 

We  now  use  the  notation  x,  y,  z  for  the  spatial  co¬ 
ordinates  and  k,  l,  m  for  the  wavenumber  coordinates. 
The  background  wind  has  components  U  in  the  x  di¬ 
rection  and  V  is  the  y  direction.  The  dispersion  relation 
is 

a>  =  khN!(kl  +  m2)m,  (32) 


Fig.  5.  The  vertical  velocity  for  the  WSK  model  predicted  by  the  present  simplified  Fourier 
method  (top)  and  by  a  numerical  simulation  (bottom).  The  contour  interval  in  both  plots  is 
0.04  ms-1,  with  positive  contours  shaded  and  the  zero  contour  omitted.  Values  range  from 
—0.37  to  0.38  ms-1  for  the  Fourier  method  and  from  —0.34  to  0.38  ms-1  for  the  numerical 
model.  The  Fourier  solution  is  calculated  at  heights  of  every  kilometer  from  the  ground  to 
an  altitude  of  30  km.  At  each  height,  256  horizontal  wavenumbers  are  used.  The  corre¬ 
sponding  spatial  domain  has  a  horizontal  periodic  length  of  about  804  km  and  a  horizontal 
grid  spacing  of  about  1.57  km.  The  numerical  model  has  130  points  in  the  horizontal  and 
52  points  in  the  vertical,  with  a  grid  spacing  of  1  km  in  both  directions. 
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Fig.  6.  The  vertical  velocity  at  z  —  2.5  km  computed  by  the  present  simplified  Fourier 
method  (top)  and  by  a  numerical  model  (bottom).  The  contour  interval  is  0.02  ms-1,  with 
positive  contours  shaded  and  the  zero  contour  omitted.  Values  range  from  —0.18  to  0.21  for 
the  present  method,  and  from  —0.21  to  0.25  for  the  numerical  model.  A  sponge  is  used  in 
the  numerical  model  beyond  about  x  =  80  km. 


where  kh  =  ( k 1  +  l2)m.  For  stationary  waves  u>  =  —kU 
-  IV. 

The  uniform  approximation  (25)  is  applicable  to  the 
three-dimensional  case  after  the  simple  modification  of 
including  the  dependence  on  /  in  h(k,  l ),  m(k,  l,  z),  etc., 
and  in  the  inverse  Fourier  transform 

r](x,  y,  z)  =  J  J  f ){k,  l,  z)eiVa+ly)  dk  dl.  (33) 

Similarly,  the  expression  (21)  for  those  Fourier  com¬ 
ponents  without  turning  points  applies  to  the  three-di¬ 
mensional  case  when  the  /  dependence  is  included. 

For  an  example  we  consider  the  WSK  model  with 
parameters  given  by  (30)  and  V  =  0.  The  topography 
is  now  three-dimensional: 

h(x,  y )  =  h0a2/(x 2  +  y2  +  a2)312,  (34) 

with  Fourier  transform 

1 

h(k,  l)  =  — h0a2e-k*a.  (35) 

27 T 

All  Fourier  components  have  turning  points,  except 
those  with  k  =  0,  for  which  u>  =  0.  We  ignore  these 
components. 

The  vertical  velocity  is  contoured  in  Fig.  6  at  a  height 
of  2.5  km.  The  top  of  Fig.  6  shows  the  solution  from 
the  present  Fourier  method,  and  the  bottom  shows  the 


results  of  a  numerical  simulation.  The  numerical  model 
is  the  same  one  used  to  produce  the  results  in  Fig.  5, 
with  the  additional  detail  that  the  grid  spacing  in  y  is  1 
km,  as  it  is  in  x  and  z.  For  the  Fourier  method,  256 
wavenumbers  are  used  in  k  and  /.  The  corresponding 
spatial  domain  has  a  horizontal  periodic  length  of  about 
670  km  and  a  horizontal  grid  spacing  of  about  1.31  km. 

6.  Final  comments 

The  usual  Fourier  integral  representation  for  moun¬ 
tain  waves  in  a  height-dependent  background  has  been 
simplified  by  replacing  the  exact  vertical  eigenfunctions 
with  their  ray  approximation.  An  Airy  function  correc¬ 
tion  removes  the  singularity  associated  with  caustics  and 
a  small  amount  of  damping  (in  the  form  of  an  imaginary 
part  of  the  wavenumber)  removes  the  singularity  as¬ 
sociated  with  resonant  modes.  The  main  result  is  given 
by  (25),  whose  inverse  Fourier  transform  is  the  spatial 
mountain  wave  solution. 

Because  of  its  reliance  on  the  ray  approximation,  the 
present  method  does  not  account  for  partial  reflection 
from  gradual  changes  in  the  background.  In  addition, 
the  Airy  function  representation  used  here  fails  in  some 
cases.  For  example,  if  the  mean  wind  in  the  model  of 
section  4  were  replaced  by  a  jet-shaped  profile,  then 
caustics  near  the  tip  of  the  jet  would  not  be  described 
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by  an  Airy  function  but  by  a  parabolic  cylinder  function 
[or  Weber  function;  see  Kravtsov  and  Orlov  (1999)]. 

The  effects  of  the  earth’s  rotation  have  also  been  omit¬ 
ted.  Shutts  (2001)  includes  rotation  in  a  mountain  wave 
model,  using  a  Fourier  method  with  the  exact  vertical 
eigenfunctions  for  a  particular  mean-wind  profile.  One 
computational  issue  is  that  an  infinitely  long  train  of 
near-inertia  waves  forms  in  the  lee  of  the  mountain, 
which  gives  periodic  wrap-around  errors  in  the  discrete 
approximation  of  the  inverse  Fourier  transform.  This  is 
the  same  problem  we  have  encountered  here  with  the 
resonant  trapped  modes  (see  the  discussion  surrounding 
Fig.  4),  and  Shutts  uses  a  similar  fix.  He  adds  a  small 
imaginary  component  to  the  intrinsic  frequency,  which 
damps  the  near-inertia  waves  and  limits  their  horizontal 
extent.  We  have  used  this  damping  and  the  ray  approx¬ 
imation  for  Shutts’  vertical  eigenfunctions  to  approxi¬ 
mately  reproduce  Fig.  6  of  Shutts  (2001),  but  we  have 
not  made  more  thorough  tests  with  rotation. 
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APPENDIX 


The  ray  solution  with  U,  and  tp, ,  corresponds  to  the 
ray  incident  on  the  caustic,  while  the  ray  solution  with 
U2  and  (I)  ,  corresponds  to  ray  reflected  from  the  caustic. 
For  our  problem,  the  incident  and  reflected  rays  have 
the  same  amplitude  dependence  on  height,  so  I  U,  I  = 

I  U2 1 ,  but  there  is  a  phase  shift  of  +  tt/2  upon  reflection 
(section  3).  Thus,  U2  =  U,  exp(/7r/2),  or  iUl  =  U2, 
which  gives  B  =  0.  For  our  problem, 

4>t  =  —  (f>2  =  —  l  \m\  dz,  (A7) 


where  we  have  written  the  phase  relative  to  the  caustic 
height  zt-  This  gives  £  =  0,  and  the  uniform  approxi¬ 
mation  reduces  to 


ii  =  -f(27T)1/2(-p)1/4[/1Ai(p).  (A8) 

If  U  represents  the  vertical  displacement,  then  by  con¬ 
servation  of  wave  action  I  U,  I  is  proportional  to  G~m, 
where  G  =  N2cgJu>  was  introduced  in  (12).  It  remains 
only  to  apply  the  boundary  conditions  at  z  =  0.  We 
could  do  that  for  the  ray  solution  U  or,  more  directly 
now,  for  the  uniform  approximation  u.  The  result  is  (25); 
that  is, 


Gjp  Ai(p) 
G2p0  Ai(p0) 


(A9) 


The  zero  subscripts  indicate  evaluation  at  z  =  0,  and  h 
is  the  Fourier  transform  of  the  topography. 


Derivation  of  the  Uniform  Approximation  (25) 

We  follow  the  method  of  Kravtsov  and  Orlov  (1999, 
section  10.1.3).  The  ray  solution  is  of  the  form 

U  =  uxe^  i  +  U2e< (Al) 

and  the  uniform  approximation  is  of  the  form 

u  =  irm[AAi(p)  +  iBAi'(p)]e‘(.  (A2) 


Here,  Ai  is  the  Airy  function  and  the  prime  means  dif¬ 
ferentiation  with  resect  to  p.  The  Airy  function  is  a  good 
approximation  of  the  linear  solution  near  the  caustic. 
The  idea  of  the  uniform  approximation  is  to  modify  the 
argument  and  amplitude  of  the  Airy  function  so  that  it 
also  fits  the  linear  solution  at  distances  farther  away 
from  the  caustic. 

A  comparison  of  (Al)  with  the  far-field  asymptotic  form 
of  (A2)  leads  Kravtsov  and  Orlov  to  the  following: 


2 

3 


-((/>!  +  (f>2). 

(A3) 

-(</>!  -  <£>-,), 

(A4) 

— p)1/4(ft/1  +  U2)e-M\ 

(A5) 

-p)m(iUl  -  Lye--'4. 

(A6) 
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